• We will use functions from packages base, utils, and stats (pre-installed and pre-loaded)
• We may also use the packages below (specifying package::function for clarity).

# Load pckgs for this R session

# --- General 
library(here)     # tools find your project's files, based on working directory
library(dplyr)    # A Grammar of Data Manipulation
library(skimr)    # Compact and Flexible Summaries of Data
library(magrittr) # A Forward-Pipe Operator for R 
library(readr)    # A Forward-Pipe Operator for R 
library(tidyr)    # Tidy Messy Data
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax

# ---Plotting & data visualization
library(ggplot2)      # Create Elegant Data Visualisations Using the Grammar of Graphics
library(ggfortify)     # Data Visualization Tools for Statistical Analysis Results
library(scatterplot3d) # 3D Scatter Plot

# --- Statistics
library(MASS)       # Support Functions and Datasets for Venables and Ripley's MASS
library(factoextra) # Extract and Visualize the Results of Multivariate Data Analyses
library(FactoMineR) # Multivariate Exploratory Data Analysis and Data Mining
library(rstatix)    # Pipe-Friendly Framework for Basic Statistical Tests
library(car)        # Companion to Applied Regression
library(ROCR)       # Visualizing the Performance of Scoring Classifiers

# --- Tidymodels (meta package)
library(rsample)    # General Resampling Infrastructure  
library(broom)      # Convert Statistical Objects into Tidy Tibbles

• Logistic regression is a type of Generalized Linear Model (GLM): a more flexible version of linear regression that can work also for categorical response variables or count data (e.g. poisson regression).

• Logistic regression is used to model a binary response variable: e.g. yes|no (or benign|malignant in the biopsy dataset).

• Logistic regression fits the coefficients \(\beta_0\), \(\beta_1\), …, \(\beta_k\) to the data using the Maximum Likelihood Estimation method (unlike linear regression, which uses the Least Squares Estimation method).

If we have predictor variables like \(x_{1,i}\), \(x_{2,i}\), …, \(x_{k,i}\) and a binary response variable \(y_i\) (where \(y_i = 0\) or \(y_i = 1\)), we need a “special” function (or "link" function) able to transform the expected value of the response variable into the \([0,1]\) outcome we’re trying to predict.

If \(p_i\) is the probability that \(y_i = 1\)… \[ logit (p_i) = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \cdots + \beta_k x_{k,i} \]

…the logit function is defined as:

\[ logit(p_i) = \ln\left( \frac{p_i}{1-p_i} \right) \]

Name: heart_data.csv
Documentation: Toy dataset prepared for teaching purposes. See reference on the data here Data Analytics with R
Sampling details: This dataset contains 10,000 observations on 4 variables.

# Use `here` in specifying all the subfolders AFTER the working directory 
heart_data <- read.csv(file = here::here("practice", "data_input", "05_datasets",
                                      "heart_data.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

Then we examine the relationship between the binary outcome variable heart_disease and the binary predictor variable coffee_drinker.

• we use the handy count function from dplyr to count occurrences in categorical variable(s) combinations.

# Dataset manipulation
heart_train %>% 
  # count the unique values per each group from 2 categorical variables' combinations
  dplyr::count(heart_disease, coffee_drinker, name = "count_by_group") %>% 
  dplyr::group_by(coffee_drinker) %>% 
  dplyr::mutate(
    total_coffee_class = sum(count_by_group),
    proportion = count_by_group / total_coffee_class) %>% 
  dplyr::ungroup() %>% 
  # filter only those with heart disease
  dplyr::filter(heart_disease == "Yes_HD") %>% 
# Plot  
  ggplot(aes(x = coffee_drinker, y = proportion, fill = coffee_drinker)) + 
  geom_bar(stat = "identity") + 
  scale_fill_manual(values = c("#57b7ff", "#e07689")) + 
  theme_minimal() +
  ylab("Percent with Heart Disease") +
  xlab("Coffee Drinker (Y/N)") +
  ggtitle("Figure 3: Percent of Coffee Drinkers with Heart Disease") +
  labs(fill = "Coffee Drinker") + 
  scale_y_continuous(labels = scales::percent)

We could technically use linear regression (Figure 1) to study likelihood of having Heart Disease in relation to risk factors, but it would not work:

• 🤔 How should we interpret a continuous prediction when it is different from \(0\) or \(1\)?

# Hypothetical linear regression model
lin_reg_model <- lm(heart_disease ~ fast_food_spend + coffee_drinker, data = heart_data)
# Plot
ggplot(heart_data) +
  geom_point(aes(fast_food_spend, heart_disease, color = factor(heart_disease)), 
             alpha = 0.25, size = 1.5) + 
  scale_color_manual(values = c("0" = "#005ca1", "1" = "#9b2339")) +
  # regression line
  geom_abline(slope = coef(lin_reg_model)[2], 
              intercept = coef(lin_reg_model)[1], color = "black", size = 1.2, linetype = "dashed" ) + 
  labs(title = "Heart disease [0,1] against fast food spend (US$) + regression line (black dashed)") + 
  ylab("Heart Disease (Y/N)") + 
  xlab("Annual Fast Food Spend") 

• Also, recall that linear regression models implies certain assumptions:
  • Linear relationship between Y and the predictors X1 and X2
    
  • Residuals must be 1) approximately normally distributed and 2) uncorrelated
  • Homoscedasticity: residuals should have constant variance
  • Non collinearity: predictor variables should not be highly correlated with each other

• Assumptions are not met in this case

# diagnostic plots
par(mfrow=c(2,2))
plot(lin_reg_model)

• Table 1 (next) shows the model output, with the coefficient estimate for each predictor.
  • The broom::tidy function converts the model summary into a more readable data frame.
  • The odds ratio (= exponentiated coefficient estimate) is more interpretable than the coefficient itself.

# Convert model's output summary into data frame
heart_model_coef <- broom::tidy(heart_model) %>% 
  # improve readability of significance levels
  dplyr::mutate('signif. lev.' = case_when(
    `p.value` < 0.001 ~ "***",
    `p.value` < 0.01 ~ "**",
    `p.value` < 0.05 ~ "*",
    TRUE ~ ""))%>%
  # add odds ratio column
  dplyr::mutate(odds_ratio = exp(estimate)) %>%
  dplyr::relocate(odds_ratio, .after = estimate) %>%
  dplyr::mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  # format as table
  knitr::kable() %>% 
  # reduce font size
  kable_styling(font_size = 20) %>% 
  # add table title
  kableExtra::add_header_above(c("Logistic Regression Analysis of Heart Disease Risk Factors"= 7))

heart_model_coef

• Intercept: gives the log-odds of heart disease when all predictor variables are zero. While not our focus, the highly negative value suggests very low baseline probability of heart disease in the sample.
• Income: Based on a p-value = 0.8182, we conclude that income is not significantly associated with heart disease. (Anyhow, the odds ratio of ≈1 would suggest no change in odds based on income.)

• Coffee Drinker (=YES): the negative coefficient of -0.7296 (log-odds) means that coffee drinking is associated with lower odds of having heart disease. Specifically, it decreases the odds ratio of heart disease by e^0.7296 which means (relative to non-coffee drinkers) an \(OR =\) 0.48 \(=\) 48%, holding all other predictors fixed. This effect is statistically significant as p-value = 0.0122.
• Fast Food Spend: the positive coefficient of 0.0024 (log-odds) means that fast food spend is associated with increased odds of having heart disease. Given this is a continuous variable, this means that for every 1-unit increase in fast food spending, the odds of heart disease increase by (e^0.0024 - 1) \(=\) (1.0024 - 1), i.e. by 0.24%. This effect is highly significant as p-value = 0.0000.

  • But we should consider a more adequate scale! (relative to fast food annual spending):

    • Δ+ = 1 USD: \((1.0024)^{1} \approx 1.0024 \quad \text{hence (1.0024 - 1)x100 = 0.24% increase in odds}\)

    • Δ+ = 10 USD: \((1.0024)^{10} \approx 1.024 \quad \text{hence (1.024 - 1)x100 = 2.4% increase in odds}\)

    • Δ+ = 100 USD: \((1.0024)^{100} \approx 1.27 \quad \text{hence (1.27 - 1)x100 = 27% increase in odds}\)

Name: Biopsy Data on Breast Cancer Patients
Documentation: See reference on the data downloaded and conditioned for R here https://cran.r-project.org/web/packages/MASS/MASS.pdf
Sampling details: This breast cancer database was obtained from the University of Wisconsin Hospitals, Madison from Dr. William H. Wolberg. He assessed biopsies of breast tumours for 699 patients up to 15 July 1992; each of nine attributes has been scored on a scale of 1 to 10, and the outcome is also known. The dataset contains the original Wisconsin breast cancer data with 699 observations on 11 variables.

The biopsy data contains 699 observations of 9 continuous variables, V1, V2, …, V9.

The dataset also contains a character variable: id, and a factor variable: class, with two levels (“benign” and “malignant”).

# check variable types
str(biopsy)

'data.frame':   699 obs. of  11 variables:
 $ ID   : chr  "1000025" "1002945" "1015425" "1016277" ...
 $ V1   : int  5 5 3 6 4 8 1 2 2 4 ...
 $ V2   : int  1 4 1 8 1 10 1 1 1 2 ...
 $ V3   : int  1 4 1 8 1 10 1 2 1 1 ...
 $ V4   : int  1 5 1 1 3 8 1 1 1 1 ...
 $ V5   : int  2 7 2 3 2 7 2 2 2 2 ...
 $ V6   : int  1 10 2 4 1 10 10 1 1 1 ...
 $ V7   : int  3 3 3 3 3 9 3 3 1 2 ...
 $ V8   : int  1 2 1 7 1 7 1 1 1 1 ...
 $ V9   : int  1 1 1 1 1 1 1 1 5 1 ...
 $ class: Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...

There is one incomplete variable V6 = “bare nuclei” with 16 missing values.

• remember the package skimr for exploring a dataframe?

# check if vars have missing values
biopsy %>% 
  # select only variables starting with "V"
  skimr::skim(starts_with("V")) %>%
  dplyr::select(skim_variable, 
                n_missing)

# A tibble: 9 × 2
  skim_variable n_missing
  <chr>             <int>
1 V1                    0
2 V2                    0
3 V3                    0
4 V4                    0
5 V5                    0
6 V6                   16
7 V7                    0
8 V8                    0
9 V9                    0

We can decide what to do in these cases (informed by our knowledge of the dataset):

• Option 1) We drop the observation with incomplete data (i.e. with missing values for V6 = “bare nuclei”) with 16 missing values.

# remove rows with missing values
biopsy_drop <- biopsy %>% 
  dplyr::filter(!is.na(V6))

mean(biopsy_drop$V6)

[1] 3.544656

• Option 2) We impute the missing values with the mean of the variable V6 = “bare nuclei”.

# impute missing values with the median of the variable
biopsy_impute <- biopsy %>% 
  dplyr::mutate(V6 = ifelse(is.na(V6), median(V6, na.rm = TRUE), V6))

mean(biopsy_impute$V6)

[1] 3.486409

• Biopsied cells of 700 breast cancer tumors, used to determine if the tumors were benign or malignant.
• This determination was based on 9 characteristics of the cells, ranked from 1(benign) to 10(malignant):
  • 1) Clump Thickness – How the cells aggregate. If monolayered they are benign and if clumped on top of each other they are malignant
  • 2) Uniform Size – All cells of the same type should be the same size.
  • 3) Uniform Shape If cells vary in cell shape they could be malignant
  • 4) Marginal Adhesion – Healthy cells have a strong ability to stick together whereas cancerous cells do not
  • 5) Single Epithelial Size – If epithelial cells are not equal in size, it could be a sign of cancer
  • 6) Bare nuclei – If the nucleus of the cell is not surrounded by cytoplasm, the cell could be malignant
  • 7) Bland Chromatin – If the chromatin’s texture is coarse the cell could be malignant
  • 8) Normal Nucleoli – In a healthy cell the nucleoli is small and hard detect via imagery. Enlarged nucleoli could be a sign of cancer
  • 9) Mitosis – cells that multiply at an uncontrollable rate could be malignant

• The explanatory variable(s) (clump_thickness, …, mitosis) can be renamed for better readability
• The observations with missing values (bare_nuclei) are removed for simplicity
  
• Patient ID (id) can be dropped as it is not used in this analysis

# Create a clean version of the dataset
biopsy_clean <- biopsy %>%
  # rename the columns (new = old)
  rename(
    id                = ID,
    clump_thickness   =  V1,
    uniform_size      =  V2,
    uniform_shape     =  V3,
    marginal_adhesion =  V4,
    single_epith_size =  V5,
    bare_nuclei       =  V6,
    bland_chromatin   =  V7,
    normal_nuclei     =  V8,
    mitosis           =  V9,
    class             =  class) %>% 
  # remove rows with missing values
  na.omit(bare_nuclei) %>% 
  # remove the id column
  select(-id)

# check the structure of the dataset
paint::paint(biopsy_clean)

data.frame [683, 10]
clump_thickness   int 5 5 3 6 4 8
uniform_size      int 1 4 1 8 1 10
uniform_shape     int 1 4 1 8 1 10
marginal_adhesion int 1 5 1 1 3 8
single_epith_size int 2 7 2 3 2 7
bare_nuclei       int 1 10 2 4 1 10
bland_chromatin   int 3 3 3 3 3 9
normal_nuclei     int 1 2 1 7 1 7
mitosis           int 1 1 1 1 1 1
class             fct benign benign benign benign benign ma~

• We fit a logistic regression model to the biopsy_train data using:
  • the glm function with argument family = binomial to specify the logistic regression model;
  • and with Class ~ . to specify an initial model that uses all the variables as predictors (backward elimination approach).

# Building initial model 
model = stats::glm(class ~ . , family = binomial, data=biopsy_train)

• Table 2 shows the model summary, with the coefficient estimate for each predictor.
  • the broom::tidy function converts the model summary into a data frame.

broom::tidy(model) %>% 
  mutate('Sign.lev' = case_when(
    `p.value` < 0.001 ~ "***",
    `p.value` < 0.01 ~ "**",
    `p.value` < 0.05 ~ "*",
    TRUE ~ ""))%>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  knitr::kable() 

In logistic regression, the coefficients are in the form of the natural logarithm of the odds of the response event happening (i.e. \(Y_i = 1\)):

\[logit(p_i) = \ln\left( \frac{p_i}{1-p_i} \right) = -9.5063 + 0.3935 \times Clump\_Thickness + ... + 0.5065 \times Mitosis\]

However, with some algebraic transformation, the logit function can be inverted to obtain the probability of the response event happening as a function of the predictors:

\[p_i = \frac{1}{1 + e^{-(-9.5063 + 0.3935 \times Clump\_Thickness + ... + 0.5065 \times Mitosis)}}\]

This equation represents the logistic regression model’s best-fit line.

• The model includes all variables, but we could make it more parsimonious by removing variables that are not significant!
• We use a statistic called the Akaike Information Criterion (AIC) to compare models.
  • AIC’s calculation gives a penalty for including additional variables.
  • The model with the lowest AIC is considered the best.

# For example let's fit a model without the variable `uniform_size`
model2 = glm(class~ .-uniform_size, family = binomial, data=biopsy_train)

According to the AIC values, the model2 seems better (AIC is lower).

# Compare the AIC values of the 2 models
tibble(Model = c("model", "model2"), 
       AIC = c(AIC(model), AIC(model2) )) %>% 
  kable()

• The MASS package’s function stepAIC enables to perform a systematic model selection (by AIC):
  • The direction argument specifies the direction of the stepwise search.
  • The trace argument (if set to TRUE) prints out all the steps.
• The best_model has removed these variables:
  • uniform_shape
  • single_epith_size
  • normal_nuclei
• The best_model has the lowest AIC value (from 100 to 98.5), despite a higher Residual Deviance than the full model (from 80 to 80.6), albeit by a very slight amount.

# Select the best model based on AIC
best_model <- MASS::stepAIC(model, direction = "both", trace = FALSE)

# Compare the AIC values of full and best model
tibble(Model = c("model", "best_model"), 
       AIC = c(AIC(model), AIC(model2)),
       Deviance = c(deviance(model), deviance(best_model))) %>% kable()

A 3D scatterplot of observations shows the first 3 principal components’ scores.

• For this one, we need the scatterplot3d() function of the scatterplot3d package;
• The color argument assigned to the Label variable;
• To add a legend, we use the legend() function and specify its coordinates via the xyz.convert() function.

# 3D scatterplot ...
plot_3d <- with(PC_scores, 
                scatterplot3d::scatterplot3d(PC_scores$PC1, 
                                             PC_scores$PC2, 
                                             PC_scores$PC3, 
                                             color = as.numeric(Label), 
                                             pch = 19, 
                                             main ="Figure 3: 3D Scatter Plot", 
                                             xlab="PC1",
                                             ylab="PC2",
                                             zlab="PC3"))

# ... + legend
legend(plot_3d$xyz.convert(0.5, 0.7, 0.5), 
       pch = 19, 
       yjust=-0.6,
       xjust=-0.9,
       legend = levels(PC_scores$Label), 
       col = seq_along(levels(PC_scores$Label)))

Biplots have two key elements: scores (the 2 axes) and loadings (the vectors). As in the scores plot, each point represents an observation projected in the space of principal components where:

• Biopsies of the same class are located closer to each other, which indicates that they have similar scores referred to the 2 main principal components;
• The loading vectors show strength and direction of association of original variables with new PC variables.

As expected from PCA, the single PC1 accounts for variance in almost all original variables, while V9 has the major projection along PC2.